First-principles calculation of mechanical properties of Si (001) nanowires and 

comparison to nanomechanical theory 
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We report the results of first-principles density functional theory calculations of the Young's 
modulus and other mechanical properties of hydrogen-passivated Si (001) nanowires. The nanowires 
are taken to have predominantly {100} surfaces, with small {110} facets according to the Wulff 
shape. The Young's modulus, the equilibrium length and the constrained residual stress of a series 
of prismatic beams of differing sizes are found to have size dependences that scale like the surface area 
to volume ratio for all but the smallest beam. The results are compared with a continuum model 
and the results of classical atomistic calculations based on an empirical potential. We attribute 
the size dependence to specific physical structures and interactions. In particular, the hydrogen 
interactions on the surface and the charge density variations within the beam are quantified and 
used both to parameterize the continuum model and to account for the discrepancies between the 
two models and the first-principles results. 

PACS numbers: 68.35.Gy, 62.25.+g, 85.85.-(-j 
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I. INTRODUCTION 

The mechanical response of structures at the nanoscale 
is known to be different than that of their macroscopic 
analogs, and surface effects in these high surface-to- 
volume devices are important^ Significant strides for- 
ward have been made in the understanding of these 
effects, but a predictive theory of nanomechanics re- 
mains an open problem both at the academic level and 
in terms of implementation in nanodevice design codes. 
Nanoscale mechanical devices have been proposed for ap- 
plications for a range of nanoelectromechanical systems 
(NEMS). These devices include high-frequency oscilla- 
tors and filters;^ nanoscale surface probes?^ probes of 
single molecules- and spins, ^ nanofluidic valves,^ and q- 
bits for quantum computation.^'^ The process of design 
and fabrication of these devices is extremely challenging, 
requiring new techniques for synthesis, characterization, 
integration and modeling of the device performance that 
are currently the subject of active research. The process 
is complicated in part by uncertainties about how even 
perfectly fabricated nanoscale mechanical devices should 
behave due to our incomplete theoretical understanding. 

At sufficiently small sizes the fact that materials are 
comprised of atoms and are not continuous media be- 
comes important. That size scale is quite close to atomic 
dimensions. Somewhat larger nanoscale structures may 
be described by continuum mechanics provided the the- 
ory is suitably extended to account for the occurrence of 
effects irrelevant in larger structures.^ One class of these 
effects is related to surfaces. Since the surface area to vol- 
ume ratio grows as the size of a structure is decreased, 
surfaces are expected to play a more prominent role at 
the nanoscale. 

Of the various mechanical properties, the Young's 
modulus is of particular interest as an important param- 
eter in the function of nanoscale devices such as flexural- 
mode mechanical resonatorsS' and as an archetype in 



terms of the scaling behavior of a variety of mechani- 
cal properties including the Poisson ratio, the anelastic 
damping coefficients and so on. The Young's modulus 
is defined as the ratio of the stress applied to stretch a 
cylindrical (or prismatic) beam to the resulting elonga- 
tion strain. For bulk materials it may be expressed in 
terms of the bulk elastic constants Cijki and it is a mate- 
rial constant, independent of the size of the structure. If 
the system is comparable in size to any mechanical inho- 
mogeneities it contains, such as grains or inclusions, the 
the modulus may exhibit a size dependence; interestingly, 
at the nanoscale even a system free from internal inho- 
mogeneities has been predicted to have a size-dependent 
Young's niodulusi^iiS 

Here we present in detail the results of an ab-initio 
study of the mechanical properties of silicon nanowires. 
Our goal is to calculate these properties from first prin- 
ciples, using a quantum-mechanical description of the 
electronic binding that is free from any empirical input 
or other a priori assumptions about the nature of the 
bonds. We then compare the results to existing nanome- 
chanical models of the size dependence of the properties 
of nanosystems. This comparison requires the calcula- 
tion of structural and mechanical properties of various 
reference systems, which we also report here. Some of 
the nanowire results have been presented previously in a 
more concise form.^^ 

While our principal focus is on the Young's modulus, 
we also calculate and analyze the residual stress and equi- 
librium elongation of the nanowires. We consider pris- 
matic Si [001] nanowires with a combination of {100} 
and {110} hydrogen-passivated surfaces, and single crys- 
tal cores as in experiment j^i^^ We have chosen the [001] 
orientation for the longitudinal axis because of its rele- 
vance to the NEMS devices;^ Si nanowires grown rather 
than etched typically have different orientations.^^ Hy- 
drogen passivation results from rinsing the oxidized Si 
surfaces with HF. It provides a standard system suitable 
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for a systematic study of the size dependence in nanome- 
chanics. With other surface conditions the band gap can 
vary greatly, and nanowires can go from semiconduct- 
ing to metallic;^* whereas the H-passivated wires remain 
semiconducting^^ down to the smallest sizes studied here 
and the surfaces do not change the nature of Si-Si chem- 
ical bonding from its covalent character. 

The article is organized as follows. We begin with 
a review in Section of the state of our current un- 
derstanding of the mechanics of nanoscale structures fo- 
cusing on nanowires, especially nanowires composed of 
Group IV elements. Then we report in Section IIIII the 
result of calculations of Si bulk and surface systems that 
are used later in the analysis of Si nanowire mechani- 
cal properties. In Section IIVI we describe the geometry 
and configuration of the nanowires studied here. We de- 
scribe how the series of wires of various sizes was created, 
and note subtleties in the description of the wires with 
continuum quantities. In Section |V] we present and ana- 
lyze the residual stress of the nanowires constrained from 
longitudinal relaxation. We analyze this stress and the 
related property of the equilibrium length of the wire. Fi- 
nally, we report the result of calculation of the Young's 
modulus, and analyze the size dependence. The basic 
calculations of the Young's modulus, the residual stress 
and the equilibrium length were reported briefly in Ref. 
[ill . Here we present the results for the reference systems 
and the detailed analysis of the mechanics not included 
in the short article. This analysis goes beyond the suc- 
cessful comparison of the first-principles results with the 
scaling laws derived from continuum mechanics models 
including surface effects presented previouslyii and al- 
lows the determination of the specific physical structures 
and interactions at the atomic and electronic level that 
lead to the size dependence. Many, but not all, of these 
structures and interactions are present in the surfaces 
of the slab reference systems that may be analyzed at 
much less computational cost than the nanowires. The 
first-principles calculations allow us to assess in detail the 
validity and point of breakdown of the continuum models 
and the physics captured therein. 



II. A BRIEF REVIEW OF NANOMECHANICS 

The mechanical properties of nanowires are expected 
to depend on the size (diameter) of the wire because 
surface effects increasingly dominate as the devices are 
miniaturized down to the nanoscale. This expectation 
has been borne out by computer simulation using atom- 
istic techniques based on empirical potentials. The first 
such calculations were done for single-crystal a-quartz 
beams, finding a systematic size dependence in which 
the Young's modulus decreased with decreasing sizei^ii^ 
These and calculations of the Young's modulus for vari- 
ous other materials have predicted a size-dependent mod- 
ulus with an additive correction to the bulk value that 
scales like the surface area to volume ratio ^i^J^ 



Continuum-based formalisms for nanoscale mechanics 
have been proposed that include the effect of surface 
properties on the mechanical behavior. One class of mod- 
els is based on the surface free energy and its first two 
strain derivatives (the surface stress and the surface elas- 
tic constants) JiiSiii^ Other approaches make use of the 
Cauchy-Born rule applied to surfaces to avoid the need 
for precalculated surface properties for metals^ and a 
Kirchoff rod model for the properties of torqued amor- 
phous nanowires Beyond the inclusion of surface prop- 
erties, there have been efforts to explore the relevance of 
other sources of size-dependence. A few studies claim 
an additional contribution that scales like the edge to 
volume ratio (cf. Ref. 0): such a contribution, with a 
factor of the logarithm of the separation of the edges, 
has been discussed for epitaxial quantum dotsi^i^i^i An- 
other study proposes the size dependence of the Young's 
modulus due to the anharmonicity (nonlinearity) of the 
bulk elastic moduli together with the strain resulting 
from the surface stress 

An intuitive way of understanding these effects is that 
there is a layer of material at the surface (and edges) 
whose mechanical properties differ from those of the bulk 
including different elastic moduli and eigenstrains. The 
term eigenstrain here means that the surface layer is con- 
strained by its interface to the bulk to be at a nonzero 
strain with respect to its minimal energy state; i.e., at 
a nonequilibrium lattice constant if the surface layer is 
crystalline. The surface layer could be chemically dis- 
tinct from the bulk, such as an oxide layer or a hydrogen- 
passivated surface, but the effect may be entirely due to 
the structural difference at the surface, such as a bare 
reconstructed surface. These surface effects could pos- 
sibly force the system to deviate significantly from the 
bulk equilibrium, and go out of the linear elastic regime. 
Atomistic calculations provide insight into these effects, 
but then we must assess to what extent the results from 
empirical potentials can be trusted. For many multi- 
component systems, empirical potentials are either not 
available or not validated. Even for pure Si nanowires, 
there is physics missing from empirical potentials that 
may be crucial, possibly to be included in the future in 
generalized potentials such as done recently for platinum 
nanoparticles.'^^ One example is the buckled dimer recon- 
struction on the Si {100} surface that is not the ground 
state of any of the standard Si empirical potentials,— or 
even a tight-binding model. 

To date, experimental data on the size dependence of 
nanostructure mechanics are very limited. If fabricating 
the nanoscale structures and measuring their mechanical 
properties such as the Young's modulus is difficult, then 
the difficulty is compounded in obtaining reproducible 
measurements free from systematic error across a series of 
structures of decreasing size. Promising work has begun 
in this direction. Atomic force microscopy (AFM) mea- 
surements of the Young's modulus^^ (E) of cast metal- 
lic nanowires with diameters in the range of 30 to 250 
nm show a strong size dependence j2£ Recent experiments 
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TABLE I; Bulk elastic properties calculated with DFT. The 
Young's modulus (E) and the Poisson's ratio (i/) have been 
derived from the elastic moduli Cu and Ci2- All units are in 
GPa except for Poisson's ratio, which is dimensionless. 
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C44 
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80.4 


Ebulk 


122.8 


125.1 
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I'bulk 


0.27 


0.26 


0.28 



have also found a strong size dependence for E of ZnO 
nanowires with diameters in the range of 17 to 550 nm^ 

For semiconductor wires, Measurements of the Young's 
modulus and the bending modulus of crystalline boron 
nanowires with diameters of 40-58 nm and 43-95 nm, re- 
spectively, have shown no systematic size dependence.- ^ 
The Young's modulus of single crystal germanium 
nanowires with the diameters of 40-160 nm is also found 
to be comparable to the bulk valuei^^ A study using 
a different AFM technique reported a value of E of 
18 ± 2 GPa for irregularly shaped < 10-nm Si (100) 
nanowires'r^ for Si(lll) wires with 100-200 nm in diam- 
eter, E has been found to be consistent with the bulk 
value*^ Si (110) nanowires with diameters of 12-170 
nm show size-dependent softeningi^ In another Group- 
IV system, measurements of E for silica nanobeams 
have demonstrated that the way in which the beam is 
clamped (i.e., the boundary conditions) affects the ap- 
parent valuci^l Experimental challenges measuring the 
intrinsic nanoscale Young's modulus make this a topic 
of continued activity, leveraging earlier work on the me- 
chanics of nanotubesj22, 

In the absence of definitive experimental data, first- 
principles quantum mechanical calculations based on 
density functional theory (DFT) can provide a robust 
prediction for the behavior of the nanoscale structures, 
but there have been few results reported. True first- 
principles techniques do not rely on any empirical data, 
solving the quantum-mechanical Kohn-Sham equations 
of DFT to achieve predictions from first principles. '^^ One 
quantum study based on an empirical tight-binding tech- 
nique has been publishedi ^^i'^^ Recently the size depen- 
dence of the Young's modulus of thin slabs has been re- 
ported calculated from first-principles (^i*^ results that 
are quite relevant but not equivalent to nanowire calcu- 
lations due to the absence of edges. We are not aware 
of any ah initio calculations of nanowire moduli in the 
literature apart from the brief article, Ref. [111. 



III. BULK AND SURFACE REFERENCE 
CALCULATIONS 

Nanowires have structural aspects ranging from bulk- 
like atomic arrangement in the core of wires to the more 



open, and often significantly relaxed, surface and edge 
structures. As the size of the wire decreases, the surfaces 
and edges play an increasing role. In this section we es- 
tablish the reference properties of the bulk silicon and the 
surfaces needed to analyze the mechanics of hydrogen- 
passivated Si (001) wires. The reference data will help 
us understand the complex mechanics of nanowires in 
terms of simpler physics and assess how continuum sur- 
face physics parameterized by the properties of the ref- 
erence systems can augment bulk continuum mechanics 
to provide a robust description of nanometer-scale struc- 
tures. 



A. Bulk properties 

A series of calculations has been done to obtain the 
bulk properties of crystalline silicon. The history of den- 
sity functional theory investigations of silicon is long, go- 
ing as far back as the late 1970's.^'^ The elastic constants 
of silicon were calculated from first principles in the work 
by Nielsen and Martin^ and anharmonic effects in sili- 
con have been mentioned in the early work by Ihm and 
coworkers. Our purpose in performing similar calcula- 
tions in the context of this study is to provide an assess- 
ment of accuracy of our results compared to literature 
values, and to provide reference numbers calculated using 
the same code and techniques as in the nanowire calcu- 
lations in order to provide directly comparable numbers 
to analyze the nanowire results. 

We use first-principles density functional theory 
(DFT) : specifically, the Vienna ab-initio simulation pack- 
age (VASP) using the projector augmented-wave method 
(PAW) *^^'^^ within the generalized gradient approxima- 
tion (GGA) by Perdew and coworkers*"'^ The PAW po- 
tentials with 4 valence electrons (3s^3p^) are used. The 
energy cutoff for the plane-wave expansion is 29.34 Ry, 
and 12x12x12 Monkhorst-Pack mesh**^ is used for k- 
point sampling. The system consists of the 8 atoms of a 
single Si diamond cubic unit cell with periodic boundary 
conditions. The supercell is deformed to the appropriate 
strain and the atomic positions are fully relaxed for each 
calculation. 

From the bulk crystal under uniaxial loading, i.e., 
stressed along the [001] direction and stress-free in two 
transverse directions, the Young's modulus and Poisson's 
ratio of the bulk reference system have been obtained. 
The equilibrium Young's modulus calculated from deriva- 
tives of a fit of the total energy in the strain range from 
5% compression to 5% tension is 122.5 GPa using a fifth- 
order polynomial fit, the corresponding Poisson's ratio 
is 0.27. The value does not change significantly for a 
lower order fit. The difference between the second order 
fit and the fifth order fit is less than 1%. The negligi- 
ble modulus difference between harmonic approximation 
and anharmonic expansion implies that the effect of bulk 
anharmonicity on the Young's modulus is small at these 
strains. Another way to assess the effect of anharmonic- 
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FIG. 1: Hydrogen passivation on Si(lOO) surface, (a) 2x1 
monohydride, (b) 1x1 symmetric dihydride, and (c) 1x1 
canted dihydride. 

ity is to examine the compression-tension asymmetry of 
the modulus using the fifth order fit: 5% tension in the 
[001] direction results in mere 3% anharmonic increase 
in the Young's modulus. It is well known that the cova- 
lently bonded silicon has a weak anharmonicity compared 
to metals, as typically expressed in terms of a relatively 
small Griineisen parameter^^ We show below based on 
these reference calculations that the anharmonicity does 
not play an important role in the size dependence of the 
Young's modulus for silicon nanowires, but the conclu- 
sion might have been different in a more strongly anhar- 
monic systemi^ 

We have also computed Cn, C12 and C44 from sepa- 
rate calculations. The values we obtain are Cn =154.6 
GPa, Ci2=58.1 GPa and C44=74.4 GPa, in good agree- 
ment with the previous DFT calculations, and are ~10% 
less than the corresponding experimental valuesi^ The 
Young's modulus and Poisson's ratio have been alterna- 
tively obtained from Cn and C12, and listed in Table [J 
The direct numbers from strictly uniaxial stress, and in- 
direct numbers from cubic elastic constants are essen- 
tially identical. 

B. Surface properties 

We now turn to the properties of the H-terminated 
Si surfaces so important to the nanowire mechan- 
ics. The ground states of low-index Si surfaces have 
long been studied both for bare (clean) reconstructed 
surfaces^i^i^ and for hydrogen-passivated surfacesi^i^ 
In recent years, after the ground states had been iden- 
tified, the attention has switched to the mechanical 
properties of the bare surfaces for the application to 
Si nanowires using empirical potentialsjii and to Si 
thin slabs using first principles calculations!^ Despite 
considerable attention to Si surfaces, the mechanics of 
hydrogen-passivated surfaces has not been studied in de- 
tail. 

The mechanics of the hydrogen-passivated surfaces, to- 
gether with the mechanics of the bare surfaces, give an 
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FIG. 2: Surface energies of the symmetric and canted H- 
terminated Si(lOO) surface and the H-terminated Si(llO) sur- 
face as a function of hydrogen chemical potential as calculated 
in DFT. The hydrogen chemical potential is calculated rela- 
tive to molecular hydrogen dissociation energy. 

indication of the range in which the surface mechanics 
in the real world resides. The two surface conditions 
represent two opposing ideals, i.e., perfect passivation of 
very reactive bonds and no passivation of such dangling 
bonds. This idealization motivates the need to under- 
stand the mechanics of the hydrogen-passivated surfaces 
qualitatively as well as quantitatively as a basis for fur- 
ther study of the mechanics of nanowires, not to mention 
many other NEMS devices. In principle, the hydrogen- 
passivated surfaces still have a surface effect, similar to 
those from the bare wires, that arises from having a free 
surface. This is an intrinsic effect. In addition, they 
have another effect due to the surface hydrogen atoms, an 
extrinsic effect. The surfactant-induced stress has been 
explained within the context of local electronic environ- 
ment for semiconductor surfaces,—"^ but it is primarily 
due to the interaction between adsorbate and substrate 
atoms. Hence, it is less important for hydrogenated 
Si surfaces, where the hydrogen- hydrogen interaction is 
prominent as evidenced in the surface (ground-state) 
structures rSS^ The efforts to understand the surface me- 
chanics using experimental®^ and simulated^ nanoinden- 
tation provide qualitative descriptions but the accurate 
(size-dependent) Young's modulus or the surface elastic 
constants have proved difficult to access via this indirect 
technique, leading to a need for direct measurement. In 
the work presented here first-principles calculations are 
crucial to deal properly with quantum mechanical effects 
such as the hydrogen-hydrogen interaction, and to inves- 
tigate directly the size-dependent mechanical properties 
of the hydrogenated surface in a quantitative sense. 

Thus, we investigate on the surface properties of the 
low-index facets relevant to (001) nanowires. Si {100} 
and {110} surfaces. A few surface patterns are known for 
the hydrogen-passivated {100} surface^ and illustrated 
in Fig. [1] we have focused on the 1x1 dihydride phase 
partially because some of the wires in the range of interest 
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TABLE II: Surface energies, stresses, and elastic constants for Si {100} surface with symmetric and canted dihydride phases, 
and Si {110} surface calculated with DPT. For surface elastic constants, [001] is taken to be the principal direction for Sn, 
and the other in-plane direction orthogonal to [001] is taken to be the second direction for 822- The surface modulus, S, is 
analogous to the bulk modulus, and defined as S — {Sn + S22 + 2S'i2)/4 when isotropic strain is applied, i.e., en = £22. 





surface energy 
(meV/A2) 


surface stress 
(meV/A2) 


^ii 


surface elastic constants (eV/A^) 

5*22 S12 


S 


{100} symmetric 


28.5 


-123.2 


-1.191 


-1.191 1.919 


0.364 


{100} canted 


19.6 


-55.0 


-0.659 


-0.659 0.457 


-0.101 


{110} 


8.2 


-1.3 


-1.223 


0.354 -0.614 


-0.526 



may be too small to develop any larger pattern fully, and 
partially because our focus in this work is on size eflFects, 
and thus we must have the same pattern for all sizes. It 
is also known that the dihydride phase is preferred over 
the monohydridc phase at low temperaturci^ 

The surface energies of hydrogenated Si surfaces we 
calculated in DFT are shown in Fig. [21 with the three 
curves corresponding to the {100} symmetric dihydride 
surface, the {100} canted dihydride surface, and the 
{110} surface. Again VASP was used, with 29.34 Ry 
cutoff energy and 12x 12x 1 Monkhorst-Pack mesh in the 
two-dimensional Brillouin zone. The calculations were 
performed on slabs with 14 layers of Si for the {100} 
slabs and 15 layers for the {110} slabs with a single unit 
cell in plane, for a total of 14 Si and 4 H atoms, and 30 Si 
and 4 H atoms, respectively. The atomic positions were 
relaxed until the force on each atom was less than 10^^ 
eV/A. 

The surface energies depend on the chemical poten- 
tial of hydrogen, and the hydrogen chemical potential 
is taken relative to the hydrogen molecule dissociation 
energy; i.e., zero chemical potential implies that the hy- 
drogens dissociate from the molecule and attach to the 
surface to terminate the dangling bonds without any cost 
in energy. The zero-point-energy terni^ is omitted, but 
this term only shifts the energy and does not alter the 
physics. The choice of the reference chemical potential is 
somewhat arbitrary but this choice measuring the chem- 
ical potential relative to the H2 dissociation energy fa- 
cilitates comparison of the two {100} surfaces: a change 
in the hydrogen chemical potential adds the same off- 
set to both surface energies, and only the surface energy 
difference matters. The horizontal line represents zero 
surface energy, and below this line the surface energy is 
negative; i.e., surface area is maximized at the cost of 
material cohesion. We are interested in the region above 
the zero-surface-energy line, and we find that the {110} 
surface is preferred over the other two surfaces. 

The surface stress can be unambiguously determined 
since it is a strain derivative of the surface energy, 
and hence the chemical potential dependence goes away. 
The surface stress has been calculated for various ma- 
terials previously using first-principles^ and classical 
atomistic^'* techniques. As can be seen in Table HIl the 
surface stress is twice larger for the symmetric {100} sur- 
face stress than for the canted {100} surface. The strong 



H-H repulsion between neighboring hydrogen atoms is 
relaxed by tilting the dihydride, and the surface stress 
is reduced. The negative stress indicates compression. 
The surface stress of the {110} surface is essentially neg- 
ligible. Across a range of loading conditions and fitting 
procedures, its values are small, varying little and what 
variation there is most likely comes from a numerical ar- 
tifact: the surface stress under [001] strain ranges from 
-1.3 to -3.2 meV/A^ depending on the order of fitting. 
The resulting uncertainty due to the fitting in the equi- 
librium surface lattice spacing is less than lO"'^ A for the 
15-layer {110} slab. 

The next derivative (second strain derivative) of the 
surface energy is the surface elastic constant. Similar to 
bulk elastic constants, surface elastic constants account 
for the material stiffness. Three surface elastic constants, 
•Sii, S22 and have been calculated with the principal 
direction of [001]. The negative constants indicate soften- 
ing due to the surface, and the opposite signs of 6*11 and 
5i2 can be thought of as a negative Poisson effect. For 
example, a higher positive 6*12 for the symmetric {100} 
surface means that, when the surface slab is stretched in 
the [001] direction, the slab would expand (direction of 
energy reduction) in the transverse direction as well due 
to the particular alignment of the dihydride in the [110] 
direction. On the other hand, a positive Poisson effect 
has been observed for the {110} slab. This conventional 
result is expected as there is virtually no H-H repulsion 
on that facet. 

Also presented is the surface modulus, 5*, a two dimen- 
sional counterpart of the bulk modulus under hydrostatic 
loading. As expected, the {100} symmetric surface has 
the highest surface modulus, the {100} canted surface is 
next, and the {110} surface shows the lowest modulus: 
the hydrogen repulsion is highest for the {100} symmet- 
ric surface, and it is virtually zero for the {110} surface. 
The strong repulsion of the symmetric dihydride even 
leads to a positive surface modulus, implying that, in the 
case of isotropic plane stress the {100} symmetric surface 
slab is even stiffer than the bulk system with the same 
volume. However, the increase in individual elastic con- 
stants shows a different tendency than that of the surface 
moduH: the 6*11 is highest for the {100} canted surface. 
We believe that this may be explained by the relaxation 
of the hydrogen repulsion during uniaxial strain. The 
uniaxial strain in the [001] direction induces the shear 
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strain for the {100} surface unit cell, and dihydrides be- 
come misaligned and consequently the associated repul- 
sion is more or less relaxed. In this situation, the H-H 
repulsion may be substantially relaxed with respect to 
the applied strain for the {100} canted surface. On the 
other hand, {100} symmetric surface seems to undergo 
a smaller change in the H-H interaction energy and its 
contribution to 5*11 is also smaller. Again, this change 
has nothing to do with the absolute magnitude of hy- 
drogen interaction. Rather, the magnitude of modulus 
is directly related to the interaction energy change with 
respect to the given strain, in particular the curvature 
with respect to the applied strain. 

To summarize, the stiffening effect of the surface hy- 
drogen for the {100} canted surface is equally noticeable 
under uniaxial and biaxial strain, but the effect seems 
prominent under biaxial strain for the {100} symmet- 
ric surface. In other words, an arbitrary strain would 
relieve the hydrogen repulsion for the {100} canted sur- 
face, but a uniform biaxial strain would dominantly do it 
for the {100} symmetric surface. This complexity of the 
hydrogenated {100} surface comes from the fact that the 
principal crystallographic directions of the bulk are (001) 
directions, whereas those of the {100} dihydride surface 
are (110) directions: the principal directions from one 
perspective are the maximum shear directions from the 
other. 

The two hydrogenated {100} surfaces, symmetric and 
canted, have substantial differences in their surface stress 
and surface moduli, and the resulting equilibrium surface 
lattice spacings for the 14-layer slab (~1.96 nm thick) 
are 1.83% and 0.76% longer than the bulk lattice, respec- 
tively. The 15-layer {110} slab (--2.95 nm thick) exhibits 
the elongation on the order of 0.01%, i.e., less than 10""^ 
A difference from the bulk spacing, due to its negligible 
intrinsic surface stress. 



C. Silane chains 

In order to separate the hydrogen interatomic (H-H) 
interaction contribution from the other contributions to 
the surface elastic constants, the effect of surface hydro- 
gen is evaluated from a periodic chain of silane molecules, 
in which the dominant H-H repulsive interaction is com- 
ing from neighboring silane molecules. The periodic 
chain is taken from the fully relaxed {100} surface slab as 
illustrated in Fig. [Sj and the backbonds are terminated 
with 2 additional hydrogen atoms. The bond angles of 
the relaxed slab geometry are preserved, and the back- 
bond length is chosen to be 1.49 A based on the structure 
of an isolated SiH4 molecule. The canted silane chain is 
shown, but the same procedure is applied to the sym- 
metric chain. The extrinsic contribution to the surface 
elastic constants from the H-H interaction is given ap- 
proximately by 




FIG. 3: (Color online) One dimensional periodic chain of 
silanes and its relationship to the (100) surface. The silane 
chain, (a), has been taken from the fully-relaxed hydrogenated 
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FIG. 4: (Color online) SiH4-SiH4 interaction energy calcu- 
lated in DFT as a function of strain for the symmetric and 
canted chains. The solid curve is an exponential fit to the 
presented data. 



where U^~^ is the H-H contribution to the total energy 
and A the area covered by one silane. The prefactor 
is due to the decomposition into the longitudinal and 
transverse directions of the wire including the Poisson 
effect with v « I'buik = 0.27. U^~^ can be approximated 
by the silane-silane interaction energy in the chain and 
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given as 

U"-"{e) « [AT^ e]/iv - (2) 

where {/^^""^[iv] is the total energy of the iV-molecule 
chain of SiH4 , and JJ'^'^^"'"'^ is the energy of a single SiH4 
molecule. U^~^ can be interpreted as a energy penalty 
to bring repulsive silane molecules into a single chain. 
It is a reasonable approximation in the sense that the 
underlying physics is essentially identical in both cases: 
H-H exchange repulsion is the dominant interaction, and 
both the SiH2 group on the {100} surface and SiH4 in the 
silane chain have almost identical local environments. 

U^^^ as a function of strain is obtained from a series 
of silane chains, whose geometries are taken from the 
corresponding relaxed {100} surface slabs. Due to the 
nature of H-H repulsion, U^~^ is fitted to an exponential 
function around the bulk lattice spacing as shown in Fig. 
m The fit of an exponential form to the data yielded the 
following dependence: 

f^symmetric = 0.359 exp(-10.97e) (3) 
C^cl"tfd = 0.167 exp(-7.10e) (4) 

in units of eV. The ideal Sii^^ obtained from the silane 
chains using Eq. [T] for the symmetric {100} and canted 
{100} slabs are 1.488 eV/A^ and 0.291 eV/A^, respec- 
tively. 

IV. GEOMETRY OF NANOWIRES 

The cross-sectional shape of the Si (001) wires we study 
is a truncated square with four {100} facets and four 
{110} facets. Some wires studied here have no {100} 
facets; for those that do, the ratio of the facet areas is 
taken to be roughly in accordance with the Wulff shape 
for a bare wire with {100}-p(2x2) and {110}-(lxl) sur- 
face reconstructions; i.e., the ratio of {100} to {110} area 
is 3.5:1. 

The reason for using the Wulff shape of bare nanowires 
is that the effect of surface condition on the material stiff- 
ness can be addressed by comparing hydrogen passivated 
wires and bare wires. In the Wulff shape of hydrogen- 
passivated wires, the {110} area would be larger than 
the {100} area due to the lower energy of that facet, op- 
posite to the Wulff shape of bare nanowires used here. 
It is advantageous to use the bare-wire shape not only 
because the surface energy of a bare slab can be unam- 
biguously determined whereas that of a slab with surface 
adsorbates such as hydrogen depends on the adsorbate 
chemical potential, but because having a common cross 
section allows us to focus on the size effect or the role 
of surface condition apart from any effect due to shape. 
The detailed comparison of the bare and H-terminated 
wires will be given elsewhercjS^ 

The electronic structure of nanowires is one dimen- 
sional. It requires a large electronic excitation to probe 




A ' Ail+ldr/r)/ 



FIG. 5: Size effect of the definition of the cross-sectional area. 
The dotted circle in the second circle indicates the initial size 
whose cross-sectional area is A. 

the transverse dimensions of the wire. The mechanics of 
the nanowire, while not fully three dimensional, has un- 
avoidable three dimensional aspects to it. The nanowires 
have a finite cross-sectional area that enters into many 
of their mechanical properties and because of this cross 
section they are able to support mechanical bending mo- 
ments. The relevance of the transverse dimensions to 
mechanical properties poses a challenge. The mechan- 
ical properties are most succinctly phrased in terms of 
continuum mechanics, but this requires the definition of 
continuum measures such as the cross-sectional area and 
the transverse dimensions in terms of atomistic proper- 
ties such as ionic positions or electron densities. 

To calculate the Young's modulus, for example, we de- 
fine the cross-sectional area to be the area bounded by 
the centers of the outermost (H) atoms. This choice is 
motivated by the fact that the volume excluded by the 
beam from access by outside atoms is determined from 
the forces arising from electron interactions. In other 
words, any experimental measure of volume is based on 
the probe being strongly repelled due to electron inter- 
actions rather than between those of nuclei. Therefore, 
a system that can be thought of as a discrete system in 
an atomic description is continuous in an electronic de- 
scription. It is unusual, however, to apply the techniques 
of continuum mechanics at subatomic length scales,— 
and indeed conventional scale-free continuum mechan- 
ics breaks down at nanometer scales or higher because of 
the lack of physics relevant to small structures such as 
surface energies and surface stresses. 

The surface definition used here ensures that most 
of the electron density is enclosed by the surface, the 
boundary formed by H atoms, and the electron den- 
sity from Si atoms essentially vanishes beyond this point. 
In addition, the positions of the nuclei are well defined 
and not subjective. Other definitions of the bounding 
surface exist: for example, the midplane between two 
identical H-passivated surfaces at their minimum energy 
separation!^ 

Had we taken a different definition of the cross- 
sectional area, the apparent size dependence of E would 
have been different, an ambiguity associated with de- 
scribing discrete atomic systems with continuum mechan- 
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ics. To illustrate the issue, consider the effect on the size 
dependence of E from changing the definition of cross- 
sectional area by modifying the position of the surface r 
by Sr as shown in Fig. O e.g. Sr would be the atomic 
radius for a surface going through atomic centers vs. one 
circumscribing their electron cloud. The change dr takes 
Ato A' ^ A{1 + 2Sr/r) 



E{r) 
E'{r) 



F 
F 



^(l + 2(5r/r)e 
£;(l-2(5r/r + ---) 



(5) 

(6) 

(7) 
(8) 



where F is the force carried by the beam and e is the re- 
sulting strain. The leading change to the Young's modu- 
lus, —2E(po) 5r/r, is proportional to the surface area to 
volume ratio. This change modifies the value coefhcient 
of the surface-area-to- volume ratio term in the Young's 
modulus, so when we calculate this coefficient, it is im- 
plicit that its value is with respect to our prescription for 
the location of the surface. The uncertainty vanishes as 
the system size is increased, i.e., the surface to volume ra- 
tio is decreased, and this level of ambiguity is completely 
negligible for the modulus of macro-scale structures. 



V. NANOWIRE CALCULATIONS 

For each of the nanowire geometries shown in Fig. [6l 
the Si atoms were initially positioned at their bulk lattice 
sites and hydrogen atoms were added to terminate the 
bonds at the surfaces, and this configuration was relaxed. 
The surface structure, such as the ground state canted 
dihydride {100} surface, is obtained naturally from the 
relaxation and requires no initial seeding. The supercell 
size of each wire is one cubic unit cell long along the wire, 
and has more than 10 A vacuum space in the transverse 
directions. The numbers of Si atoms and H atoms in the 
supercell for each of the geometries are also shown in Fig. 

m 

Again the first-principles DFT code VASP has been 
used for the calculations, as described above. The en- 
ergy cutoff for the plane- wave expansion is 29.34 Ry, and 
1x1x12 Monkhorst-Pack meshjiS or six points in the one 
dimensional irreducible Brillouin zone, is used for fc-point 
sampling. The relaxation convergence criterion was that 
the force on each atom be smaller than 2 x 10~'^ eV/A for 
the 1.49-nm and smaller wires, and 10^2 gV/A for the 
larger wires, for which convergence of the Young's mod- 
ulus was attained with the less stringent tolerance. The 
residual stress, equilibrium length and Young's modulus 
have been calculated for the range of wire sizes shown in 
Fig.H 

The extent of the range of wire sizes is dictated by 
computer resource limitations. The calculations for the 



largest wire, the 3.92-nm wire, required roughly 922 
hours using typically 128 processors on the 2.3 GHz 
Xeon-based MCR supercomputer, for a total of 109 285 
CPU-hours for this one wire. For the final stages of the 
relaxation, in which the calculation setup as described 
above was applied, it was necessary to use up to 512 pro- 
cessors. 



A. Residual stress 

Residual stress can be a problem or a feature of doubly 
clamped oscillators. During etching in the device fabri- 
cation process, stress forms in the mechanical structure. 
This stress shifts the resonant frequency of the mechani- 
cal beam, and it can affect other properties such as dissi- 
pation. For example, recently resonators with very high 
quality factors have been demonstrated at room temper- 
ature using silicon nitride nanobeams under high tensile 
stressi^ Residual stress may arise due to a number of 
physical phenomena. One example is the residual stress 
coming from surface stress due to interactions in the pas- 
sivation layer, and it is this effect in the H-passivated Si 
nanowires that we investigate now. Within our DFT cal- 
culations, the residual stress is equated to the axial stress 
when the longitudinal wire lattice spacing is fixed at the 
bulk value, (Jzz{.Lq). 

To calculate the residual stress, the system was ini- 
tially relaxed to its zero-temperature minimum energy 
with the length of the periodic supercell held fixed at the 
bulk lattice spacing in the longitudinal direction. Then 
the relaxed total energy was calculated for each beam in 
a series of longitudinal strains at increments of roughly 
0.5%. These total energy values were fit to a polynomial. 

We then obtain the residual stress and other important 
properties from the strain derivatives of the polynomial: 
the first derivative gives the axial stress, Uzz, the min- 
imum of the polynomial, Czz-eq, gives the strain corre- 
sponding to the equilibrium length, and the value of the 
curvature at the minimum gives the Young's modulus E: 



CTzzi^zz) 
^zz {^zz — eq^ 



V- 





^dU/dez 



E = V-^d^U/de 



(9) 
(10) 

(11) 



where U is the DFT total energy. The equilibrium length 
Leg is determined from the equilibrium strain and the 
initial length Lq according to 

Leg = io (1 + e 



zz~eq ) 



(12) 



where here and throughout the article we are using en- 
gineering strain, which is adequate for the small strains 
of interest. The details of the equilibrium length and the 
Young's modulus are discussed in the next two subsec- 
tions. We focus on the stress in this subsection. 

The size dependence of the residual stress of the H 
passivated Si nanowires evident in Fig. [7] is driven by 
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405Si lOOH 




3.92 nm 

FIG. 6: (Color online) Cross sections of fully relaxed hydrogen-passivated wires, with each Si atoms colored as shown in the 
legend corresponding to its transverse relaxation in A. The numbers above the wires stand for the number of atoms in the 
supercell. For example, 405Si lOOH means that the supercell has 405 Si atoms and 100 H atoms. The numbers below represent 
the wire width, where the width is defined as the square root of the cross-sectional area. 



compressive surface stress. The residual stress may be 
decomposed into core and surface contributions. The 
latter may be further decomposed into extrinsic surface 
contributions from hydrogen (H-H) interactions and in- 
trinsic surface contributions from the change to the Si 
bonds near the surface compared to the Si bulk (Si-H 
and modified bond order Si-Si). Since DFT only pro- 
vides a total energy, this decomposition is somewhat am- 
biguous. Techniques such as Bond Order Potentials^ 
and the Generalized Pseudopotential Theory^^ have been 
introduced as formalisms for extracting atomic interac- 
tions from the underlying quantum mechanics. They are 
not sufficiently developed to apply to nanowires, however. 
Local moment techniques also provide a partition of the 
total energy atom by atom, but they too are insufhciently 
accurate for our purposes. Instead, we separate different 
contributions to the nanowire properties through the use 
of reference systems. In the first instance, the continuum 
models require core and surface properties: the uniformly 
strained bulk Si systems of Section fill Al provide an ap- 
proximation to the core of the nanowire and the slab 
surface calculations of Section IIII Bl give an approxima- 
tion to the nanowire surfaces. In a more refined analysis 
we decompose the surface contribution further and esti- 
mate the extrinsic surface (H-H) interactions to be equal 
to those of neighboring hydrogens in two silane molecules 
in the orientation and separation of the H-passivated sur- 
face. The intrinsic contribution is the remainder, i.e., the 
part that does not come from the core or the extrinsic 
surface interactions. The further decomposition of the 
surface effects into intrinsic and extrinsic enables us to 
begin to attribute the size dependence to specific physi- 
cal processes of bonding and deformation in the surface, 
sub-surface and core regions of the wire. 

The extrinsic contribution to the residual stress is the 
most important, as we now show. The intrinsic surface 
stress is small, as expected since the dangling Si bonds 
are well terminated with H atoms and the Si-Si bond or- 
der is not significantly different than in the bulk. The 



small magnitude of the intrinsic stress is best seen in 
the case of the 1.39-nm wire for which the elongation is 
less than 0.1% compared to ^--^1.5% of the 1.49-nm wire. 
The absence of {100} facets on this wire leads to a small 
extrinsic stress since the H-H separation on the {110} 
facets is relatively large. This smallness was already con- 
firmed in the previous section where the surface stress of 
the {110} surface has been found to be negligible. The 
difference between these two facets is that the vacant Si 
sites above the facets are filled by one and two H atoms on 
{110} and {100}, respectively, and the double occupancy, 
albeit with ^2 A H-H separation due to the shorter Si-H 
bond, leads to more repulsion for {100} as discussed in 
Sec. H. 

The extrinsic surface stress due to the H-H repulsion 
on the {100} facets quantitatively accounts for both the 
compressive residual stress CFzziLo) and the elongated 
equilibrium length L^q of the nanowires. They are re- 
lated to leading order through the linear elasticity: 

(^zz{Lo) = azz{core) + ^^T^'Jwi (13) 

i 

{Ln-L,q)/L,q ~ azz{Lo)/E (14) 

where A is the cross-sectional area, Wi is the width, ri^ 
is the longitudinal surface stress of facet i, and Lq is the 
bulk length of the beam. E is the Young's modulus of 
the beam. In principle Eq. (|14p has anharmonic correc- 
tions, but the strains are small and the harmonic ap- 
proximation should be good. For constant surface stress, 
the second term in Eq. (|13p is proportional to the sur- 
face area to volume ratio; the core stress is too, since 
the surface stress causes a transverse expansion of the 
wire that induces a tensile core stress. It may be seen 
in Fig. [5] that surface Si atoms on the {100} facets un- 
dergo a substantial transverse expansion, where for ex- 
ample a 0.3 A expansion amounts to a 13% change in the 
bulk Si-Si bond length. It also induces the deformation 
noticeable throughout the wire, extended deep into the 
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FIG. 7: (Color online) Silicon nanowire axial stress as a func- 
tion of wire size calculated in DFT. The virial stresses are 
obtained directly from virial formula at the given cutoff en- 
ergy indicated in the parenthesis, whereas the DFT stress is 
deduced from Eq.[9] The predictions of Eq.[l3]are also plotted 
for which the surfaces stresses are given in Table Ull 
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FIG. 8: Silicon nanowire equilibrium elongation strain as a 
function of wire size calculated in DFT. The solid curve is a fit 
to C /w of the elongation strain to four data points from 1.49- 
nm and larger wires, with C=1.9%-nm. The dotted curve is a 
fit to C /w of the elongation strain to the 1.39-nm wire, with 
C=0.14%-nm. 



core except for some of the high symmetry lines where 
the expansions toward opposite directions cancel. Hence, 
the transverse displacement would be hardly seen when 
the {100} facets are not present: the {110}-facet wires 
with no {100} facets, compared to those with the {100} 
facets, have a negligible smface expansion and the result- 
ing in-plane deformation of the core. 

The core stress arises as a kind of Laplace pressure and 
may be estimated through a generalized Young-Laplace 
equation relating the compression of the core to the sur- 
face stress. The derivation of the expression is somewhat 
lengthy and presented elsewhere, The result for the ax- 
ial stress in the core averaged across the cross section of 
the wire is 

(T,,(core) = --i.r^i"°>/w;, (15) 

TT 

where w is the width of the nanowire, v = Ci2/(Cii + 
C12) is the Poisson ratio and Cn and C12 are the bulk 
elastic constants. Here a compressive surface stress (r < 
0) leads to to a tensile axial stress {a^z > 0). 

The residual stress values evaluated with a few differ- 
ent approaches are plotted in Fig. [71 The curves denoted 
as 'virial' are the stresses directly obtained from the cal- 
culations via virial theorem. The stress curve denoted 
simply as 'DFT stress' is due to Eq. ((9]), and the one 
denoted as 'prediction' is due to Eq. The virial 

stress evaluation is poor even for reasonably high cutoff 
energies, i.e., the predicted stress in some cases has the 
opposite sign with 29.34 Ry cutoff, and a reasonable eval- 
uation requires a prohibitively high cutoff energy. On the 
other hand, the stress obtained by the strain derivative of 
the total energy with the same cutoff energy is less prone 
to the basis set incompleteness critical to the direct mea- 
sure, and shows better agreement with the stress with 



the highest cutoff energy tested. In addition, it is advan- 
tageous as we can obtain the stress information of larger 
wires, where the cutoff energy is restricted due to compu- 
tational cost. Considering the discrepancy between the 
stress due to Eq. ^ and the virial stress with the highest 
cutoff energy, the true stress is likely to be between the 
two: the convergence error for direct evaluation seems to 
be negative and polynomial fitting error may be positive. 

For the prediction based on the bulk and surface en- 
ergies, using the values from Tables [11 and Ull in Eq. ([T51) 
gives predictions in very good agreement with the full 
nanowire calculations as shown in Fig. [71 The scatter for 
1.49- and 2.05-nm wires may be accounted for by small 
edge effects. The 0.61-, 1.00- and 1.39-nm wires have no 
{100} facets and almost no residual stress as described 
above. In the case of the second smallest (0.92-nm) wire, 
all of the {100} atoms undergo substantial relaxation, as 
shown in Fig. [51 lowering the magnitude of the surface 
stress and the elongation. This high level of agreement 
gives us confidence that we understand the physics of the 
size dependence of the residual stress. 



B. Equilibrium length 

The residual stress is a direct measure of the effect of 
surface stress on doubly clamped nanowires. The anal- 
ogous property of unconstrained (free floating or can- 
tilever) nanowires is the equilibrium length. It is a prop- 
erty that is, at least in principle, directly measurable 
through x-ray diffraction. 

The equilibrium elongation plotted in Fig. [H shows a 
systematic increase in the elongation as the size of the 
wire is reduced, with the data falling into two series. The 
series are distinguished by the amount of {100} surface 
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discussed below. A solid curve representing the 
function C/w has been superimposed on one series of 
data, where w is the width of the wire. The good fit of 
this function, with C = 1.9%— nm, demonstrates that the 
elongation is proportional to the surface-area-to-volume 
ratio. The axial stress shown in Fig. [7] exhibits a similar 
trend, with the stress increasing as the size of the wire 
decreases. Again, the deviation of the 0.92-nm wire from 
the curve can be explained by the reduction in the surface 
repulsion due to the substantial relaxation around the 
edges. Equivalently, it can be explained that the wire is 
too small to have well defined {100} facets that would 
otherwise give comparable surface effects to those on the 
larger wires. 

As is the case with the residual stress (Fig. [T]), the 
wires in the series with elongation values closer to zero 
have no {100} facets. These wires show a significantly 
smaller amount of elongation compared to those in the 
series with {100} facets. The elongation for the 1.39-nm 
wire is less than 0.1%, and to achieve a comparably small 
amount of elongation in a wire with {100} facets, the wire 
would need to be 19 nm in width based on the above 
fit. The elongation would also be negligible for a larger 
wire with the same cross-sectional shape, as confirmed 
by considering that the elongation for the 2.95-nm thick 
{110} slab is less than 10"^ A. 

While the elongation of the {110}-faceted wires is con- 
siderably smaller than that of their {100}-faceted coun- 
terparts, it is actually larger than expected: as can be 
seen from the dotted curve in Fig. [SI the elongation of 
the smallest wire is even larger by roughly a factor of 2 
than the C/w extension from the 1.39-nm wire, let alone 
any potential size-dependent effects beyond the surface 
effect on the 1.39-nm wire. We do not have a conclusive 
explanation for this enhanced size dependence. It could 
be due to an edge effect, but it seems more likely that it 
is a nonlocal surface-surface interaction effect: for these 
smallest wires, each facet begins to interact with oppo- 
site facet via electronic interactions. For example, the 
1.00-nm wire has only 5 Si layers across its cross section. 
It is known from slab calculations that the effective sur- 
face properties are modified if the slab is too thin, and 
we expect an analogous effect here. 



C. Young's modulus 

We now consider the size dependence of the Young's 
modulus, the principal subject of this work. The Young's 
modulus is the second derivative of the DFT total energy 
with respect to the applied longitudinal strain divided by 
the volume pT|) . As a second derivative, it is the most 
sensitive of the quantities we compute, determining the 
values of the planewave cutoff, number of k points and 
residual force tolerance quoted above. Also, it is poten- 
tially susceptible to an unwanted dependence on the poly- 
nomial fitting procedure used for the energy. In principle, 
higher order fits should be more accurate, accommodat- 



ing any anharmonicity but at the cost of the need for 
additional data (more relaxed total energy calculations). 
The Young's modulus of the 1.49-nm wire gives an indi- 
cation of the sensitivity to the order of the polynomial fit, 
as described in Ref. 11. Fortunately, for the given order 
of fit, a higher cutoff energy does not significantly im- 
prove the fit; i.e., the largest convergence error is 0.63% 
for the 2nd-order fit compared to the highest order fit. 
We find that the second order fit with 29.34 Ry energy 
cutoff is reasonably good, differing by less than 2% com- 
pared with all the higher-order combinations tested, and 
it is the second order fit that we use for the analysis of 
all of the nanowires, and it permits direct comparison 
of the results from the entire range of nanowire sizes up 
to 405 Si and 100 H atoms at a tolerable computational 
cost. We have used this technique for a systematic study 
of the size dependence of the Young's modulus. 

The main result is that the Young's modulus becomes 
softer monotonically as the size is decreased as shown 
in Fig. [5] and that decrease is well described by a con- 
tinuum model. It drops from the bulk value (iS'JJJJ = 
122.5 GPa) roughly in proportion to the surface area to 
volume ratio. It does not exhibit a strong dependence 
on the ratio of the {100} to {110} area seen in the equi- 
librium length. The smallest Young's modulus we find 
for any of the nanowires studied here is 29.4 GPa, the 
value for the smallest (0.61-nm) nanowire. This value is 
larger than the 18 ±2 GPa reported from experiments on 
an irregularly shaped < 10-nm Si (100) nanowire.^'* Since 
the shapes are different, direct comparison is challenging, 
but we note that the experimental stress-strain curve is 
complex. Some regions of the stress-strain curve are in 
better agreement with the range of values we report here 
than the region from which the value of 18 GPa was ex- 
tracted. In any case, they do report a softening of the 
Young's modulus as we have found in the first-principles 
calculations. 

From continuum mechanics neglecting edge and non- 
local effects, the modulus can be expressed, slightly gen- 
eralizing Ref. [13, as 

i 

where S"^*^ is the surface elastic constant, a strain- 
derivative of the surface stress including both extrinsic 
and intrinsic parts. Thus, the Young's modulus may be 
decomposed into core and intrinsic, and extrinsic surface 
contributions as was done for the residual stress. 

We found above that the extrinsic H-H interactions 
dominated the surface stress and the residual stress of the 
nanowires. It is less clear a priori what should dominate 
the Young's modulus, but the fact that the modulus is 
insensitive to the facet ratio (whether the {100} facets 
are present or not) suggests several conclusions: 

• The core anharmonicity is irrelevant since the mod- 
ulus is not correlated with the equilibrium elonga- 
tion; 
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• The extrinsic contribution to the modulus (which 
is strongly facet dependent) is small; 

• The intrinsic surface elastic constant dominates 
and its {100} value may be nearly sufficient to de- 
termine E. 

Additional evidence supports each of these conclusions, 
as we now show. 

To quantify the core contribution, we calculated that 
the Young's modulus of the bulk crystal increases by only 
1.6% when strained ^^1.5% to match the most strained 
(0.92-nm) wire. This change is negligible compared to 
the observed softening (contrary to the finding that the 
bulk anharmonicity dominates the size dependence of 
the Young's modulus of embedded-atom-method copper 
nanowires^^) . 

The extrinsic effect is also small, but not negligi- 
ble. Based on silane interaction forces for the canted 
{100} surface geometry, the extrinsic contribution to the 
Young's modulus can be described as 

^E^-^'-jT.S^i'^-^i"""' (17) 

i 

where wj^^^^ is the width of a {100} facet i. We have es- 
timated that the extrinsic contribution is ^^8 GPa for the 
1.49-nm wire, roughly equal to E(1.49nm)-E(1.39nm), 
i.e., the difference in the moduli with and without {100} 
facets. Here, the canted dihydride is particularly impor- 
tant in that the initial (unrelaxed) symmetric dihydrides 
on the {100} facets of nanowires relax to a configuration 
very similar to the canted dihydrides as illustrated in Fig. 
ini only somewhat strained due to the symmetry breaking 
at the edges that destroys the in-plane periodicity of the 
infinite surface. Also, given the agreement for the resid- 
ual stress and elongation between our continuum-based 
prediction with the canted {100} surface properties and 
the DFT result, we can again confirm that the canted sur- 
face is indeed relevant to the relaxed wires. Otherwise, 
the predicted residual stress and the elongation based on 
the symmetric {100} facet would be twice larger than the 
actual value for the 1.49-nm or larger wires.^- 

The intrinsic contribution accounts for most of the 
Young's modulus. By definition, it is the remainder 
once the core and extrinsic contributions have been sub- 
tracted, and those contributions are small as we just 
showed. We may go further and create a qualitative 
map of the intrinsic contribution using a bond-strength 
calculation akin to an Einstein model with independent 
harmonic oscillators. A small longitudinal displacement 
is applied to each atom, and by measuring the induced 
force, the spring constant for each atom is deduced. Each 
atom has three spatial degrees of freedom and hence three 
oscillators, and only the longitudinal oscillators are con- 
sidered here. As can be seen in Fig. [TOl surface atoms 
have substantially softer bonds. The true intrinsic ef- 
fect might be even greater since the force built up on the 
neighboring atoms by displacing a single atom is relieved 



somewhat by the relaxation of its neighbors. The fluctu- 
ation of the spring constant for fully coordinated atoms 
may be explained by the electron density variation due to 
the surface relaxation and the natural structure of silane 
chains, i.e., Friedel-like oscillations. Some recent work 
has considered the effect of electron density variation on 
mechanical properties through a bond-order approach. — 
We have also calculated the size dependence of the 
modulus using Eq. (fTC)) based on the surface elastic con- 
stants 5''^^°*'^ from the surface reference system presented 
in Table |TT1 The results, shown in Fig. [HI are in good 
agreement with the full first-principles calculation, and 
adding the core contribution slightly improves the agree- 
ment. The reason for the scatter for 1.49- and 2.05-nm 
wires is two-fold. First, there is a variation in the tilting 
angle of the hydrides on the wire surface: those around 
the facet center are equal or closer to the symmetric con- 
figuration, i.e., no tilting or smaller tilting angle. The 
modulus estimation based only on the canted surface 
elastic constants could overshoot the real value. The size 
of the symmetric dihydride region at the middle of the 
facet is determined by the usual kink analysis and is es- 
sentially independent of the size of the facet. Second, the 
scatter for 1.49- and 2.05-nm wires may be due in part 
to small edge effects. Both the effect of the symmetric 
dihydride region and that of the edges would scale as 
the ratio of their constant areas to the facet area, and 
would therefore be less important for larger wires. Also 
plotted in Fig. [5] is the best fit curve of Ref. 17 from 
Stillinger- Weber (SW) empirical molecular statics calcu- 
lations. The SW bulk Young's modulus is 13% lower and 
the coefficient C of the 1/w term is 29% lower, repre- 
senting a weaker surface effect than in DFT. The errors 
compensate for each other, leading to reasonable agree- 
ment for the nanoscale wires. This level of agreement 
is unexpected since the SW potential does not have the 
relevant nanophysics in its functional form or its fitting 
database and the strength of the bonds does not change 
at the surface. Also, the SW nanowire calculation does 
not include a H-terminated surface, and thus the residual 
stress and equilibrium elongation are quite different. 



VI. CONCLUSIONS 

In conclusion we have found that first-principles cal- 
culation of several mechanical properties of silicon wires 
predicts a size dependence at the nanoscale. The form 
of the size dependence is in good agreement with previ- 
ous models based on empirical atomistic and continuum 
techniques, for all but the smallest wires in which effects 
such as the electronic interaction between surfaces are not 
captured in the models. The calculations presented here 
enabled an analysis of the magnitude of surface and edge 
effects in the nanowire Young's modulus from first prin- 
ciples. In each case the size dependence scales roughly 
as the surface area to volume ratio, but for different rea- 
sons. For the equilibrium length and residual stress it 
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FIG. 9: Silicon nanowire Young's modulus as a function of 
wire size as calculated in DPT. For comparison the predic- 
tions of continuum formula (|16|) are also plotted, for which the 
surface elastic constants are obtained from the {100} canted 
slabs and the {110} slabs. See Section IlIII for details. The 
solid curve E = E^Jik — C/w, with C=66.11 GPa/nm, is a 
best fit to a pure surface area to volume size dependence. 




FIG. 10; (Color online) Spring constant of Si atoms for the 
1.49-nm wire. The inset illustrates the location of each atom 
in the cross section: a quarter of the cross section is shown 
with the rest of the wire related by mirror symmetry, where 
atom number 1 is located at the center of the cross section. 

is due to the extrinsic surface stress from interactions 
in the H passivation layer; for the Young's modulus it 
arises from the intrinsic contribution to the surface elastic 
constant. The equiUbrium length and residual stress de- 
pend strongly on whether {100} facets are present or not, 
whereas the Young's modulus was essentially insensitive 
to the facet type. Surface parameters from slab calcula- 
tions capture most, but not all, of the physics. The size 
effect is not strong for the H-terminated surfaces studied 
here: the Young's modulus is softened by about 50% for 
a 1-nm diameter wire. It may be possible to measure 
this effect directly using either AFM deflection or res- 
onant frequency measurements in a double-clamped or 



cantilever configuration. The change in the equilibrium 
length is measurable with x-ray or electron diffraction 
techniques. These effect could be substantially stronger 
in silicon nanowires with different surfaces, such as bare 
or oxide surfaces, making measurement easier. These sys- 
tems are more challenging for first-principles calculation 
due to a greater number of candidate structures and a 
greater role for charge transfer in the mechanics. It is 
not clear whether the Young's modulus would increase 
or decrease as the size of the beam is reduced. There is 
much to be learned still. 

In the wires we have studied the surface atoms are pas- 
sivated by hydrogen atoms so that the chemical bonding 
is the same as in the bulk and quantum confinement is ev- 
ident, the local electronic environment of each Si atom is 
uniform throughout the wire apart from small variations. 
The lack of any noticeable quantum mechanical difference 
between the Si atoms in the core and those at the surface 
(the Si-H interface) makes it easier to apply continuum 
modeling successfully to the hydrogen-passivated wires. 
As we consider how these effects transfer to other kinds 
of wires, the extrinsic surface effect coming from the sur- 
face adsorbates will require a careful treatment, and in 
cases where the surface interactions are stronger it may 
prove more difficult for continuum models to provide an 
accurate description of small wires. 

The calculations presented here are at absolute zero 
temperature. Given the present-day computers, it is not 
possible to carry out these first-principles calculations 
at finite temperature. Based on the results of empirical 
molecular dynamics , the general form of the size 
dependence of the Young's modulus is expected to be re- 
tained at finite temperatures well below the melt point, 
but the value of the modulus will change. Naturally, 
thermal softening will shift the entire curve, but also the 
value of the coefficient of the surface-area-to- volume ratio 
term will be somewhat temperature dependent. For this 
and many other applications it is desirable to construct 
classical interatomic potentials, whether quantum-based 
or strictly empirical, that capture the surface physics 
relevant to the mechanics of nanostructuresi^ First- 
principles calculations, such as those presented here, lay 
the groundwork for the development of those potentials. 
There is clearly much to be done in the development of 
Nanomechanics . 
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